library(raster)

library(plyr)
library(dplyr)

library(ggplot2)
library(viridisLite)
library(colorspace)
library(RColorBrewer)

library(sf)
library(stringr)

memory.limit(size = 160000)
## [1] Inf

#display.brewer.pal(n = 11, name ="RdYlBGn)
colorvec <- brewer.pal(n = 11, name ="RdYlGn")

1 Data

land_mask <- raster("Masks/land_sea_mask_1degree.nc4") 
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land


elev <- raster("Worldclim/wc2.1_10m/wc2.1_10m_elev.tif") 
elev.df <- elev %>% projectRaster(to = land_mask) %>% as.data.frame(xy = T) %>% setNames(c("lon","lat", "z"))


ipcc_regions <- shapefile("Masks/IPCC-WGI-reference-regions-v4.shp") %>% spTransform(crs("EPSG:4326")) 
ipcc_regions.raster <- ipcc_regions %>% rasterize(land_mask)
ipcc_regions.df <- as.data.frame(ipcc_regions.raster, xy = T) %>% setNames(c("lon","lat","Continent","Type","Name","Acronym"))
tas.all_annual <- read.table("CMIP6/tas/tas.all_annual.txt")
pr.all_annual <- read.table("CMIP6/pr/pr.all_annual.txt")
sfcWind.all_annual <- read.table("CMIP6/sfcWind/sfcWind.all_annual.txt")
hfls.all_annual <- read.table("CMIP6/hfls/hfls.all_annual.txt")
hfss.all_annual <- read.table("CMIP6/hfss/hfss.all_annual.txt")
hurs.all_annual <- read.table("CMIP6/hurs/hurs.all_annual.txt")

cmip6_annual <- tas.all_annual %>% merge(pr.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_annual, "cmip6_annual.txt")
tas.all_january <- read.table("CMIP6/tas/tas.all_january.txt")
pr.all_january <- read.table("CMIP6/pr/pr.all_january.txt")
sfcWind.all_january <- read.table("CMIP6/sfcWind/sfcWind.all_january.txt")
hfls.all_january <- read.table("CMIP6/hfls/hfls.all_january.txt")
hfss.all_january <- read.table("CMIP6/hfss/hfss.all_january.txt")
hurs.all_january <- read.table("CMIP6/hurs/hurs.all_january.txt")

cmip6_january <- tas.all_january %>% merge(pr.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_january, "cmip6_january.txt")
tas.all_february <- read.table("CMIP6/tas/tas.all_february.txt")
pr.all_february <- read.table("CMIP6/pr/pr.all_february.txt")
sfcWind.all_february <- read.table("CMIP6/sfcWind/sfcWind.all_february.txt")
hfls.all_february <- read.table("CMIP6/hfls/hfls.all_february.txt")
hfss.all_february <- read.table("CMIP6/hfss/hfss.all_february.txt")
hurs.all_february <- read.table("CMIP6/hurs/hurs.all_february.txt")

cmip6_february <- tas.all_february %>% merge(pr.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_february, "cmip6_february.txt")
tas.all_march <- read.table("CMIP6/tas/tas.all_march.txt")
pr.all_march <- read.table("CMIP6/pr/pr.all_march.txt")
sfcWind.all_march <- read.table("CMIP6/sfcWind/sfcWind.all_march.txt")
hfls.all_march <- read.table("CMIP6/hfls/hfls.all_march.txt")
hfss.all_march <- read.table("CMIP6/hfss/hfss.all_march.txt")
hurs.all_march <- read.table("CMIP6/hurs/hurs.all_march.txt")

cmip6_march <- tas.all_march %>% merge(pr.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_march, "cmip6_march.txt")
tas.all_april <- read.table("CMIP6/tas/tas.all_april.txt")
pr.all_april <- read.table("CMIP6/pr/pr.all_april.txt")
sfcWind.all_april <- read.table("CMIP6/sfcWind/sfcWind.all_april.txt")
hfls.all_april <- read.table("CMIP6/hfls/hfls.all_april.txt")
hfss.all_april <- read.table("CMIP6/hfss/hfss.all_april.txt")
hurs.all_april <- read.table("CMIP6/hurs/hurs.all_april.txt")

cmip6_april <- tas.all_april %>% merge(pr.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_april, "cmip6_april.txt")
tas.all_may <- read.table("CMIP6/tas/tas.all_may.txt")
pr.all_may <- read.table("CMIP6/pr/pr.all_may.txt")
sfcWind.all_may <- read.table("CMIP6/sfcWind/sfcWind.all_may.txt")
hfls.all_may <- read.table("CMIP6/hfls/hfls.all_may.txt")
hfss.all_may <- read.table("CMIP6/hfss/hfss.all_may.txt")
hurs.all_may <- read.table("CMIP6/hurs/hurs.all_may.txt")

cmip6_may <- tas.all_may %>% merge(pr.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_may, "cmip6_may.txt")
tas.all_june <- read.table("CMIP6/tas/tas.all_june.txt")
pr.all_june <- read.table("CMIP6/pr/pr.all_june.txt")
sfcWind.all_june <- read.table("CMIP6/sfcWind/sfcWind.all_june.txt")
hfls.all_june <- read.table("CMIP6/hfls/hfls.all_june.txt")
hfss.all_june <- read.table("CMIP6/hfss/hfss.all_june.txt")
hurs.all_june <- read.table("CMIP6/hurs/hurs.all_june.txt")

cmip6_june <- tas.all_june %>% merge(pr.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_june, "cmip6_june.txt")
tas.all_july <- read.table("CMIP6/tas/tas.all_july.txt")
pr.all_july <- read.table("CMIP6/pr/pr.all_july.txt")
sfcWind.all_july <- read.table("CMIP6/sfcWind/sfcWind.all_july.txt")
hfls.all_july <- read.table("CMIP6/hfls/hfls.all_july.txt")
hfss.all_july <- read.table("CMIP6/hfss/hfss.all_july.txt")
hurs.all_july <- read.table("CMIP6/hurs/hurs.all_july.txt")

cmip6_july <- tas.all_july %>% merge(pr.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))
 
write.table(cmip6_july, "cmip6_july.txt")
tas.all_august <- read.table("CMIP6/tas/tas.all_august.txt")
pr.all_august <- read.table("CMIP6/pr/pr.all_august.txt")
sfcWind.all_august <- read.table("CMIP6/sfcWind/sfcWind.all_august.txt")
hfls.all_august <- read.table("CMIP6/hfls/hfls.all_august.txt")
hfss.all_august <- read.table("CMIP6/hfss/hfss.all_august.txt")
hurs.all_august <- read.table("CMIP6/hurs/hurs.all_august.txt")

cmip6_august <- tas.all_august %>% merge(pr.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_august, "cmip6_august.txt")
tas.all_september <- read.table("CMIP6/tas/tas.all_september.txt")
pr.all_september <- read.table("CMIP6/pr/pr.all_september.txt")
sfcWind.all_september <- read.table("CMIP6/sfcWind/sfcWind.all_september.txt")
hfls.all_september <- read.table("CMIP6/hfls/hfls.all_september.txt")
hfss.all_september <- read.table("CMIP6/hfss/hfss.all_september.txt")
hurs.all_september <- read.table("CMIP6/hurs/hurs.all_september.txt")

cmip6_september <- tas.all_september %>% merge(pr.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_september, "cmip6_september.txt")
tas.all_october <- read.table("CMIP6/tas/tas.all_october.txt")
pr.all_october <- read.table("CMIP6/pr/pr.all_october.txt")
sfcWind.all_october <- read.table("CMIP6/sfcWind/sfcWind.all_october.txt")
hfls.all_october <- read.table("CMIP6/hfls/hfls.all_october.txt")
hfss.all_october <- read.table("CMIP6/hfss/hfss.all_october.txt")
hurs.all_october <- read.table("CMIP6/hurs/hurs.all_october.txt")

cmip6_october <- tas.all_october %>% merge(pr.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_october, "cmip6_october.txt")
tas.all_november <- read.table("CMIP6/tas/tas.all_november.txt")
pr.all_november <- read.table("CMIP6/pr/pr.all_november.txt")
sfcWind.all_november <- read.table("CMIP6/sfcWind/sfcWind.all_november.txt")
hfls.all_november <- read.table("CMIP6/hfls/hfls.all_november.txt")
hfss.all_november <- read.table("CMIP6/hfss/hfss.all_november.txt")
hurs.all_november <- read.table("CMIP6/hurs/hurs.all_november.txt")

cmip6_november <- tas.all_november %>% merge(pr.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_november, "cmip6_november.txt")
tas.all_december <- read.table("CMIP6/tas/tas.all_december.txt")
pr.all_december <- read.table("CMIP6/pr/pr.all_december.txt")
sfcWind.all_december <- read.table("CMIP6/sfcWind/sfcWind.all_december.txt")
hfls.all_december <- read.table("CMIP6/hfls/hfls.all_december.txt")
hfss.all_december <- read.table("CMIP6/hfss/hfss.all_december.txt")
hurs.all_december <- read.table("CMIP6/hurs/hurs.all_december.txt")

cmip6_december <- tas.all_december %>% merge(pr.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_december, "cmip6_december.txt")

2 Aridity Index


cmip6 <- read.table("cmip6_annual.txt")

breaks.martonne <- c(0,5,10,20,"30","40",Inf)
cat.martonne <- c("[0,5]" = "Desert", "(5,10]"= "Arid", "(10,20]" = "Semi-arid","(20,30]" = "Temperate", "(30,40]" = "Humid", "(40,Inf]" = "Forest")

breaks.unesco <- c(-Inf,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(-Inf,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid")

#display.brewer.pal(n = 11, name ="RdYlBu")
colorvec <- brewer.pal(n = 11, name ="RdYlBu")

col.martonne <- c("Desert" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Temperate" = colorvec[8], "Humid"= colorvec[10],"Forest"= colorvec[11])
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10])

2.1 De Martonne


cmip6$AIm <- with(cmip6,pr*60*60*24*365/(tas-273.15+10)) # pr in mm/y, tas in Ceslsius


cmip6$cat.AIm <- cmip6$AIm %>% cut(breaks.martonne, include.lowest = T) %>% revalue(cat.martonne)

ggplot(subset(cmip6, period == "1970_2000")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AIm))+
  scale_fill_manual(values = col.martonne)+
  theme_void()+
  theme(legend.position = "bottom")

2.2 FAO

ea from Rh:

ea = 6.108(RH/100)exp(17.27*T/(T+237.3)) Refs https://www.nature.com/articles/s41598-023-40499-6 : Allen and Tetens

2.2.1 ET0 by month

id_vars <- c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")

cmip6_annual <- mutate(read.table("cmip6_annual.txt"),
                  t = tas - 273.15,
                  Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                  P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                  es = 0.61078*exp(17.27*t/(t+237.3)),
                  ea = (hurs/100)*es,
                  delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                  gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                  ET0_annual = ((0.408*delta*Rn+gamma*900/((t)+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind)))) 

cmip6_01 <- mutate(read.table("cmip6_january.txt"),
                   t = tas - 273.15,
                   pr_01 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                  delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                  gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                  ET0_01 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind)))) # Evapotranspiration, G considered to be negligible. in mm/day. Wind speed at 10 m is converted to wind speed at 2 m by multiplying by 0.748

cmip6_02 <- mutate(read.table("cmip6_february.txt"),
                   t = tas - 273.15,
                   pr_02 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                  delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                  gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                  ET0_02 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_03 <- mutate(read.table("cmip6_march.txt"),
                   t = tas - 273.15,
                   pr_03 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                  delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                  gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                  ET0_03 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_04 <- mutate(read.table("cmip6_april.txt"),
                   t = tas - 273.15,
                   pr_04 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                  es = 0.61078*exp(17.27*t/(t+237.3)),
                  ea = (hurs/100)*es,
                  delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                  gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                  ET0_04 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_05 <- mutate(read.table("cmip6_may.txt"),
                   t = tas - 273.15,
                   pr_05 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                  delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                  gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                  ET0_05 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_06 <- mutate(read.table("cmip6_june.txt"),
                   t = tas - 273.15,
                   pr_06 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                  delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                  gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                  ET0_06 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_07 <- mutate(read.table("cmip6_july.txt"),
                   t = tas - 273.15,
                   pr_07 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                  delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                  gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                  ET0_07 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_08 <- mutate(read.table("cmip6_august.txt"),
                   t = tas - 273.15,
                   pr_08 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                  delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                  gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                  ET0_08 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_09 <- mutate(read.table("cmip6_september.txt"),
                   t = tas - 273.15,
                   pr_09 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                  delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                  gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                  ET0_09 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_10 <- mutate(read.table("cmip6_october.txt"),
                   t = tas - 273.15,
                   pr_10 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                  delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                  gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                  ET0_10 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_11 <- mutate(read.table("cmip6_november.txt"),
                   t = tas - 273.15,
                   pr_11 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                  delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                  gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                  ET0_11 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_12 <- mutate(read.table("cmip6_december.txt"),
                   t = tas - 273.15,
                   pr_12 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                  delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                  gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                  ET0_12 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

et0_monthly <- select(cmip6_annual, id_vars) %>%
  merge(select(cmip6_01, c(id_vars, "ET0_01", "pr_01")), by = id_vars) %>%
  merge(select(cmip6_02, c(id_vars, "ET0_02", "pr_02")), by = id_vars) %>%
  merge(select(cmip6_03, c(id_vars, "ET0_03", "pr_03")), by = id_vars) %>%
  merge(select(cmip6_04, c(id_vars, "ET0_04", "pr_04")), by = id_vars) %>%
  merge(select(cmip6_05, c(id_vars, "ET0_05", "pr_05")), by = id_vars) %>%
  merge(select(cmip6_06, c(id_vars, "ET0_06", "pr_06")), by = id_vars) %>%
  merge(select(cmip6_07, c(id_vars, "ET0_07", "pr_07")), by = id_vars) %>%
  merge(select(cmip6_08, c(id_vars, "ET0_08", "pr_08")), by = id_vars) %>%
  merge(select(cmip6_09, c(id_vars, "ET0_09", "pr_09")), by = id_vars) %>%
  merge(select(cmip6_10, c(id_vars, "ET0_10", "pr_10")), by = id_vars) %>%   
  merge(select(cmip6_11, c(id_vars, "ET0_11", "pr_11")), by = id_vars) %>%
  merge(select(cmip6_12, c(id_vars, "ET0_12", "pr_12")), by = id_vars)

et0_monthly$spr <- with(et0_monthly, pr_01*31 + pr_02 * 28 + pr_03 * 31 + pr_04 * 30 + pr_05 * 31 + pr_06 * 30 +
                           pr_07 * 31 + pr_08 * 31 + pr_09 * 30 + pr_10 * 31 + pr_11 * 30 + pr_12 * 31)
et0_monthly$sET0 <- with(et0_monthly, ET0_01*31 + ET0_02 * 28 + ET0_03 * 31 + ET0_04 * 30 + ET0_05 * 31 + ET0_06 * 30 +
                           ET0_07 * 31 + ET0_08 * 31 + ET0_09 * 30 + ET0_10 * 31 + ET0_11 * 30 + ET0_12 * 31)

cmip6 <- merge(cmip6_annual, et0_monthly, by = id_vars)
cmip6$AI <- with(cmip6, abs(spr/sET0))  
cmip6$AI_annual <- with(cmip6, abs((pr*60*60*24*365)/(ET0_annual*365)))

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")

cmip6$cat.AI <- cmip6$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6$cat.AI[which(cmip6$sET0 < 400)] <- "Cold"
cmip6$cat.AI %>% unique()

cmip6$cat.AI <- cmip6$AI %>% cut(breaks = breaks.unesco) %>% 
  revalue(cat.unesco)

cmip6$cat.AI_annual <- cmip6$AI_annual %>% cut(breaks = breaks.unesco) %>% 
  revalue(cat.unesco)

View(subset(select(cmip6, c("lon","lat","sET0","AI","cat.AI","ET0_annual","ET0_annual","cat.AI_annual")), cat.AI_annual == "Hyperarid"))

for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "1970_2000" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
  scale_fill_manual(values = col.cat)+
  labs(title = i)+
  theme_void()+
  theme(legend.position = "bottom")
print(g)
}

for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "1970_2000" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI_annual))+
  scale_fill_manual(values = col.cat)+
  labs(title = i)+
  theme_void()+
  theme(legend.position = "bottom")
print(g)
}

write.table(cmip6, "cmip6_AI.txt")

cmip6 <- read.table("cmip6_AI.txt")

cmip6_AI <- cmip6 %>% subset(period == "1970_2000") %>%
   group_by(lon,lat,Continent, Type, Name, Acronym, lm, model, z) %>% 
   dplyr::summarise(tas = mean(tas, na.rm = T) - 273.15, pr = mean(pr, na.rm = T) * 60 * 60 * 24 * 365,
                   sfcWind = mean(sfcWind, na.rm = T), Rn = mean(Rn, na.rm = T),
                   ea = mean(ea, na.rm = T), ET0 = mean(sET0, na.rm = T), AI = mean(AI, na.rm = T)) %>%
  ungroup() %>%  
  mutate(source = "CMIP6")%>%
   select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn", "ea","z","ET0","AI"))


write.table(cmip6_AI, "cmip6ref.txt")

3 Maps ETO by model

cmip6 <- read.table("cmip6_AI.txt")
for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "1850_1880" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
  scale_fill_manual(values = col.cat)+
  labs(title = paste(i, "1850-1880", sep = ", "))+
  theme_void()+
  theme(legend.position = "bottom")
print(g)
}

for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "1985_2015" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
  scale_fill_manual(values = col.cat)+
  labs(title = paste(i, "1985-2015", sep = ", "))+
  theme_void()+
  theme(legend.position = "bottom")
print(g)
}

for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP245" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
  scale_fill_manual(values = col.cat)+
  labs(title = paste(i, "SSP 2-4.5, 2030-2060", sep = ", "), fill = "")+
  theme_void()+
  theme(legend.position = "bottom")
print(g)
}

for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP245" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
  scale_fill_manual(values = col.cat)+
  labs(title = paste(i, "SSP 2-4.5, 2070-2100", sep = ", "), fill = "")+
  theme_void()+
  theme(legend.position = "bottom")
print(g)
}

for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP370" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
  scale_fill_manual(values = col.cat)+
  labs(title = paste(i, "SSP 3-7.0, 2030-2060", sep = ", "), fill = "")+
  theme_void()+
  theme(legend.position = "bottom")
print(g)
}

for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP370" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
  scale_fill_manual(values = col.cat)+
  labs(title = paste(i, "SSP 3-7.0, 2070-2100", sep = ", "), fill = "")+
  theme_void()+
  theme(legend.position = "bottom")
print(g)
}

for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP585" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
  scale_fill_manual(values = col.cat)+
  labs(title = paste(i, "SSP 5-8.5, 2030-2060", sep = ", "), fill = "")+
  theme_void()+
  theme(legend.position = "bottom")
print(g)
}

for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP585" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
  scale_fill_manual(values = col.cat)+
  labs(title = paste(i, "SSP 5-8.5, 2070-2100", sep = ", "), fill = "")+
  theme_void()+
  theme(legend.position = "bottom")
print(g)
}

4 Anomalies

cmip6s <- read.table("cmip6_AI.txt") %>% group_by(lon,lat,Continent, Type, Name, Acronym, lm, model, z, period) %>% 
   dplyr::summarise(t = mean(t, na.rm = T), tas.sd = sd(t, na.rm = T), # t in °C
                    pr = mean(spr, na.rm = T) , pr.sd = sd(spr, na.rm = T),  # pr in mm/y
                   sfcWind = mean(sfcWind, na.rm = T), sfcWind.sd = sd(sfcWind, na.rm = T),
                   Rn = mean(Rn, na.rm = T), Rn.sd = sd(Rn, na.rm = T), # in MJ/m2/day
                   ea = mean(ea, na.rm = T), ea.sd = sd(ea, na.rm = T), # in kPa
                   ET0 = mean(sET0, na.rm = T), ET0.sd = sd(sET0, na.rm = T), # in mm/y
                   AI = mean(AI, na.rm = T), AI.sd = sd(AI, na.rm = T))

cmip6s$AI2 <- with(cmip6s, pr/ET0)
breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")

cmip6s$cat.AI <- cmip6s$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6s$cat.AI[which(cmip6s$ET0 < 400)] <- "Cold"

cmip6s <- read.table("cmip6s.txt")
cmip6s$anomalies <- NA
cmip6s$change.cat <- NA
cmip6s$cat.AI.ref <- NA
cmip6s$change <- NA

for(i in 1:nrow(cmip6s)){
  AIi = cmip6s$AI[i]
  lati = cmip6s$lat[i]
  loni = cmip6s$lon[i]
  AIref = subset(cmip6s, lat == lati & lon == loni & period =="1970_2000") %>% pull(AI) 
  cmip6s$anomalies[i] <- AIi-AIref
  
  catAIi = cmip6s$cat.AI[i]
  catAIrefi = subset(cmip6s, lat == lati & lon == loni & period =="1970_2000") %>% pull(cat.AI) 
  cmip6s$change.cat[i] <- ifelse(catAIi == catAIrefi, 0, 1)
  cmip6s$cat.AI.ref[i] <- catAIrefi
  cmip6s$change[i] <- paste(catAIrefi, catAIi, collapse = " ", sep = " to ")
}

write.table(cmip6s, "cmip6s.txt")
cmip6s <- read.table("cmip6s.txt")

cmip6s$dryer <- with(cmip6s, ifelse(change == "Arid to Humid", "no",
                                ifelse(change == "Arid to Semi-arid", "no",
                                ifelse(change == " Arid to Dry subhumid", "no", 
                                ifelse(change == "Arid to Hyperarid", "yes",
                                ifelse(change == "Arid to Cold", "no",       
                                ifelse(change == "Cold to Humid", "no",
                                ifelse(change == "Cold to Hyperarid", "yes",
                                ifelse(change == "Dry subhumid to Cold", "no",
                                ifelse(change == "Dry subhumid to Humid", "no",
                                ifelse(change == "Dry subhumid to Arid", "yes",
                                ifelse(change == "Dry subhumid to Semi-arid", "yes",
                                ifelse(change == "Humid to arid", "yes",
                                ifelse(change == "Humid to Semi-arid", "yes",
                                ifelse(change == "Humid to Cold","no",
                                ifelse(change == "Humid to Dry subhumid","yes",
                                ifelse(change == "Semi-arid to Dry subhumid", "no",
                                ifelse(change == "Semi-arid to Arid", "yes",
                                ifelse(change == "Semi-arid to Humid", "no",
                                ifelse(change == "Semi-arid to Cold", "no",
                                ifelse(change == "Hyperarid to Arid", "no",
                                ifelse(change == "Hyperarid to Semi-arid", "no",
                                ifelse(change == "Hyperarid to Dry subhumid", "no",
                                ifelse(change == "Hyperarid to Humid", "no",
                                NA))))))))))))))))))))))))

4.1 Difference between future and reference AI


ggplot(subset(cmip6s, anomalies < 1 & anomalies > -1))+geom_boxplot(aes(x=period, y = anomalies))


# for(i in unique(cmip6s$period)){
# g <- ggplot(data = subset(cmip6s, period == i)) + geom_raster(aes(x=lon, y = lat,  fill = anomalies))+
#   labs(title = paste("AI changes between",i, "and 1970-2000", sep = " "))+
#   theme_void()+
#   theme(legend.position = "bottom")
# print(g)
# }
for(i in c("1850_1880","1985_2015")){
g <- ggplot(data = subset(cmip6s, period == i & model == "historical" & change.cat == 1)) + geom_raster(aes(x=lon, y = lat,  fill = dryer))+
  borders()+
  scale_fill_manual(values =  c(colorvec[7], colorvec[4]), na.translate = F)+
  labs(title = paste("Is",i, "dryer than 1970-2000", sep = " "))+
  theme_void()+ylim(-55,90)+
  theme(legend.position = "bottom")
print(g)
}


for(i in c("2030_2060","2070_2100")){
  for (j in c("SSP245","SSP370","SSP585")){
g <- ggplot(data = subset(cmip6s, period == i & model == j)) + geom_raster(aes(x=lon, y = lat,  fill = dryer))+
  borders()+
  scale_fill_manual(values =c(colorvec[7], colorvec[4]), na.translate = F)+
  labs(title = paste("AI changes between",i, "and 1970-2000,", j, sep = " "))+
  theme_void()+ylim(-55,90)+
  theme(legend.position = "bottom")
print(g)
  }
}

knitr::knit_exit()